//===========================================================================
/*!
 *  \file PotHF.h
 */
//===========================================================================
#ifndef _PotHF_H
#define _PotHF_H


#include "Potential.h"
#include "comm/matrix_macro.h"

#define NPARAM 40

//#include <gsl/gsl_linalg.h>
//void LinearSolver_LU(int N,double A_data[],double x0[],double b_data[]);

//define mememo style used in class, it can be 1 for dynamic, 2 for static
#define MEM_STYLE 2



/*!
 *  \class PotHF
 *  \brief Potential models for Water Cluster
 *  \pre Need extern TTM functions to be linked
 */
class PotHF:public Potential{	
	public:
		
		PotHF(){ //!< constructor			
			LoadingDefaultParam();
		}
				
		PotHF( int dim_){//!< constructor
			#if MEM_STYLE == 1
 				r=r2=NULL;	
 				Scd=dScd=Sdd=dSdd=NULL;
 				T=D=Dtmp=Dinv=NULL;
 				E=mu=q=NULL;
 				dr=NULL;
			#endif
				
			model=HF_OSS2;
			unit=HARTREE;

			SetDim_(dim_);	LoadingDefaultParam();
		}
		
		~PotHF();
		
		void SetDim_(int dim_);		
		
		void SetExtra_(void *p){ nO=*((int*)p); //cout<<" nF = "<<nO<<endl;
		}
		
		const char* classname(){ return "PotHF";};
		
		double eval_gradient_(double const* _p, double* _g){return 0;}		
		
		int ReadParam(const char *filename); //!<read parameters for potential model from file

		int ReadParam(const double *alpha_);

		void PrintInfo(){ 
			cout<<" nAtom = "<<nAtom<<" nF= "<<nO<<" dim = "<<dim<<" equation= "<<equation()<<" unit= "<<GetUnit_()<<endl;
			cout<<" Using following param : "<<endl;
			for(int i=0;i<NPARAM;i++) cout<<alpha[i]<<endl;
		};  //!< print info for checking

		double  evaluate_( double const *_p);
		void gradient_(double const* _p, double* _g);	
		
		bool HasAnalyticalGradients() { return true; }
				
	private:
		void LoadingDefaultParam(){ //!< Loading the default parameter, done after initilization
			double t[NPARAM]={0.96000,	1.95230,	0.61680,	-0.21334,	332.2852922, 	2.84683483,	1.0956864,	0.00000204,	0.01180,	2.14617,	1.00277,	3.73356,	0.11178,	41.19416,	 0.74860,	-41.01776,	0.74859,	27.87439,	0.75515,	-1.05883,	-0.02831,	0.15970,	-0.02512,	-0.28188,	-0.36116,	-0.10125,	0.07372,	0.30578,	1.10914,	-0.11205,	-0.31964,	-0.41244,	0.56539,	-0.04931,	-1.93150,	-0.26473,	3.31874,	4.10023,	4.56354,	0.83000};

			for(int i=0;i<NPARAM;i++) alpha[i]=t[i];

			model=HF_OSS2;
			unit=HARTREE;
		};
		
	private: 
		int nAtom,nO;
		
		double alpha[NPARAM];
		double xtemp[SIZE_MAX][DGR];
		double forcetemp[SIZE_MAX][DGR];	

		double huy_oss2_energy_improve(double x[][DGR], int nAtom, int nO, double alpha[],int inverse=0);
		void huy_oss2_energy_grad(const double x[][DGR], int nAtom, int nO, double alpha[], double& energy, double grad[][3]);

		void ReleaseVars();
		
		double oss2energyRed(const double x[][DGR],int nAtom,int nO,double *alpha);
		void oss2gradient(const double x[][DGR],int nAtom,int nO,double *alpha,double grad[][DGR]);

	private: //for Huy's OSS2
	#if MEM_STYLE == 1 
		// 2 dim 	
		double **r,**r2; // distance matrices
		double **Scd,**dScd,**Sdd,**dSdd;	
		double **T,**D,**Dtmp,**Dinv;

		//1 dim 
		double *E,*mu,*q; // E,mu, charge vector
		//3 dim
		double*** dr;
		void inv(double **A, double **Ainv, int N);			
		int IR_Gauss_Seidel(int n,double **A_data,double x[], const double b[],int max_iter, double &tol);
	
	#elif  MEM_STYLE == 2
		// 2 dim 	
		double r[SIZE_MAX][SIZE_MAX],r2[SIZE_MAX][SIZE_MAX]; // distance matrices
		double Scd[SIZE_MAX][SIZE_MAX],dScd[SIZE_MAX][SIZE_MAX],Sdd[SIZE_MAX][SIZE_MAX],dSdd[SIZE_MAX][SIZE_MAX];	
		double T[SIZE_MAX][SIZE_MAX],D[SIZE_MAX][SIZE_MAX],Dtmp[SIZE_MAX][SIZE_MAX],Dinv[SIZE_MAX][SIZE_MAX];

		//1 dim 
		double E[SIZE_MAX],mu[SIZE_MAX],q[SIZE_MAX]; // E,mu, charge vector
		//3 dim
		double dr[SIZE_MAX][SIZE_MAX][3];
		void inv(double A[][SIZE_MAX], double Ainv[][SIZE_MAX], int N);			
		int IR_Gauss_Seidel(int n,double A_data[][SIZE_MAX],double x[], const double b[],int max_iter, double &tol);	
		
	#endif
		
};;

#endif